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It is shown by example that the predicted temperature dependence of the elastic constants is a useful 
measure of the ability of an effective pair potential to estimate the high temperature thermal properties of a 
metal. Our example is based on a model pair potential constructed for aluminum. This potential predicts the low 
temperature elastic constants and phonon dispersion relations with good accuracy (± a few percent). The high 
temperature elastic constants for this model potential are determined using the Monte Carlo method and are 
found to be approximately independent of temperature. Since the elastic constants of aluminum are strongly 
decreasing functions of temperature, this potential is seen to be a poor one for determining the properties of 
aluminum. We conclude that the temperature dependence of the elastic constants is a useful further test of pair 
potentials which satisfy the low temperature tests currently employed. 
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1. Introduction 

Pseudopotential theory has made it possible to develop 
effective, volume dependent pair potentials for metals [l] 1 
and thereby to estimate the high temperature thermal 
properties of metals using the methods of statistical me- 
chanics [2-4]. When a rriodel pseudopotential is used, it is 
necessary to fit the model parameters to various physical 
properties such as the low temperature lattice parameter 
and elastic constants or dispersion curves in the harmonic 
approximation. [5], Fitting an effective potential to low 
temperature properties is no guarantee that the high tem- 
perature properties of the solid will also be satisfactorily 
predicted. Cohen et al. [4] have suggested that the ability of 
an effective potential to predict reliably the high tempera- 
ture elastic constants of a metal is a good test of the poten- 
tial. In this paper we show by example that this is the case 
and assert that the predicted temperature dependence of 
the elastic constants is a useful measure of the ability of an 
effective potential to estimate high temperature thermal 
properties of a metal. 

Our discussion is in terms of a model potential for alumi- 
num which produces good dispersion curves at low temper- 
atures but which produces totally unacceptable estimates 
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for the temperature dependence of the elastic constants at 
high temperatures. The elastic constants are determined at 
high temperatures, where anharmonic effects are impor- 
tant, using Monte Carlo techniques [2-4] to evaluate the 
fluctuation expressions of Squire, et al. [2]. While the static 
lattice estimates for elastic constants employ only a discrete 
set of values of the derivatives of the potential, the full 
microscopic expressions [2] sample in principle the entire 
range of the potential when anharmonic effects ars impor- 
tant. Thus the predicted temperature variation of the elastic 
constants for a fixed volume depends on the overall shape of 
the potential. This is demonstrated by our example. 

This paper is organized as follows. Section 2 contains a 
discussion of the pair potential we have devised for alumi- 
num. We show that it satisfies the compressibility sum rule 
and state that it reliably predicts the low temperature elastic 
constants and dispersion curves for aluminum. Section 3 
describes the Monte Carlo calculation of the temperature 
dependence of the elastic constants at high temperatures 
and Section 4 summarizes our results. 



2. Pair potential 

Numerous effective pair potentials have been proposed 
for aluminum and the shape of the potential has been found 
to be quite sensitive the details of the electron screening 
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function. Duesbery and Taylor have shown that a necessary 
condition for a potential is that the screening function 
satisfy the electron gas compressiblity theorem [6], In the 
same vein, Hafner and Schmuck [7] have stated that this 
compressibility sum rule should be at least approximately 
satisfied if the correct relation between longitudinal and 
transverse branches is to be realized. 

We have found empirically that if the combination of 
elastic constants 2C' — C = C u — C l2 ~ C 44 is correctly 
predicted then the dispersion curves for the tranverse 
modes in the (1,1,0) direction are correctly ordered. Here 
Cii, C i2 and C 44 are cubic elastic constants in the usual 
notation. [8] We examined the model pseudopotential of 
Dagens, et al [9] which satisfies the compressibility sum 
rule and found that it predicts a crossing of the dispersion 
curves for the transverse modes in the (1,1,0) direction. We 
have also examined the first principles pseudopotential of 
Harrison [l]and found that it too predicts mode crossing. 
Hafner and Schmuck [7] developed a corrected Harrison 
potential and used it to calculate the elastic constants of 
aluminum at zero temperature. Their results are in agree- 
ment with the experimental values of Kamm and Ahlers 
[10]. 

The model potential used in this calculation approxi- 
mately satisfies the compressibility sum rule and cor- 
rectly predicts the value of the combination of elastic 
constants 2C' — C = C xx — C x2 ~ C 44 at zero temperature. 
Since errors which do not influence the compressibility rule 
may affect 2C" — C, we observe that a potential which satis- 
fies both conditions at zero temperature is more likely to 
function adequately at non zero temperatures than one 
which does not. As we shall demonstrate, this is not a suffi- 
cient condition. From pseudopotential theory we know that 
given the energy wave number characteristic, F(q), we can 
obtain by Fourier transformation the indirect (electron gas 
induced) interaction between the ion cores. The full poten- 
tial is the sum of the couloumb and indirect interactions. 
For simplicity we have used a model F(q) as recently applied 
to the calculation of the structure factor of liquid aluminum 
[11]. This model F(q) introduced by Ashcroft and Langreth 
[12] is expressed as 



where kj- is the Fermi wave vector for a free electron at the 
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where Ze is the charge on the ion, Q is the volume occupied 
by an ion, q is the wave vector, R lon is a parameter given a 
value of 1.12 atomic units for aluminum and e(q) is the 
dielectric function given by [13] 



density in question. The function R(q) is defined by 
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and uj(q) is the exchange screening function. We have used 
the Hubbard model for xp(q); 



w(q) = 



1 



q 2 Hk/ 



with the parameter assigned the value | = 1.699 so that the 
compressibility sum rule is satisfied. The full potential V(r) 
is given by 
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This potential is displayed in figure 1 as a function of r. It 
adequately reproduces the low temperature phonon disper- 
sion relations in the harmonic approximation [14] and also 
avoids the crossing of the nearly degenerate small k 
transverse modes in the (1,1,0) direction. Thus it satisfies 
what we consider to be a minimum requirement for a 
"good" potential function. The local maximum located bet- 
ween the 1st and 2nd neighbor distances is a common 
feature of many effective potentials for aluminum [6,9]. 

The zero temperature (static lattice) values of 2C '— C and 
the lattice compressibility /? T can be simply expressed as 
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FIGURE 1. The effective pair potential for aluminum at a density of 
(Z.Z) 2. 718 g/cm 3 (lattice parameters - .404 nm). 

The energy is in units of electron volts and the distance is in nanometers. The arrow indicates the 
place where the potential was set to zero in the Monte Carlo Calculations. 
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with r n the distance to the rc th neighbor shell and a is the lat- 
tice parameter. These expressions are obtained by summing 
over shells Squire et al's [2] equations for the elastic con- 
stants of the static FCC lattice. The values of the first five 
/,'s for the potential of figure 1 are listed in table 1. 

Table I 
n f n — dynes/cm 



20.5790 

2.0479 

-2.0199 

0.0872 

0.0515 



The predicted low temperature value of 2C — C is ob- 
tained by using table 1 and eq (2.5) to obtain 

2C -C=2.31X 10 n dynes/cm 2 . 

This is to be compared with the experimental value [15] of 
2.327 x 10 11 dynes/cm 2 . The calculated value for the com- 
pressibility is 

3//3 T = 1.694 x 10 12 dynes/cm 2 

which is to be compared with the experimental value of 
1.612 x 10 12 dynes/cm 2 . The lattice compressibility sum rule 
is satisfied to within 5 percent and the 2C'—C rule to within 
1 percent by this potential. 



3. Monte Carlo calculations 

The Monte Carlo method of Metropolis, et al [3] was used 
to calculate the thermal properties of the crystalline model 
using the effective pair potential of figure 1. The system 
consisted of 256 particles located initially on the FCC sites 



of a 4 x 4 x 4 block of unit cells with lattice parameter a = 
0.404 nm. Periodic boundary conditions were used to mimic 
an infinite system. A simulation started with an "ageing" 
period of 3X10 5 trials. The data were extracted from the 
lxlO 6 trials following the ageing sequence. The program ad- 
justed the maximum step size allowed in a trial move so that 
approximately J /2 of the proposed moves were accepted. 
Double precision arithmetic was used to accumulate the 
data and all other calculations employed single precision 
arithmetic. The use of double precision is necessary if 10 6 
samples are to be averaged using 32-bit words. 

The temperature dependent part of the pressure was 
evaluated by averaging the virial. The elastic constants were 
obtained by evaluating the expressions developed by 
Squire, et al. [2] and our program was checked by reproduc- 
ing their results for the Lennard-Jones solid. 

Estimates were obtained for the energy per particle (E), 
the pressure (p), the isothermal elastic constants (CiJ, C X2 T , 
C 4 4 7 '), the difference between the isothermal and adiabatic 
elastic constants (AC), the specific heat at constant volume 
(C v ) and the derivative of the pressure with respect to tem- 
perature at constant volume ((dP/d T)v). Fluctuation ex- 
pressions for C v and (3P/dT}v were used to estimate these 
quantities [4-16]. 

The full expressions for the pressure, C M and C i2 as func- 
tions of volume and temperature contain electron gas terms 
which are functions of the volume, but are independent of 
configuration and temperature [4]. Since we are concerned 
only with the temperature variation of these quantities, we 
have not included these contributions in our results which 
are listed in table 2. The quoted uncertainties are based on 
the stability of the estimates as the Monte Carlo process 
progresses. This is illustrated in figure 2 where the esti- 
mates for P, (dP/dT)v and C X J at 500 K are shown as a 
function of the number of trials after the ageing process had 
been completed. The program provides these estimates 
every 5X10 4 trials. Only 8.5 x 10 5 trials were made for 
r=700K due to a computer malfunction which made con- 
tinuation of that run impossible. 

Table 2 



Temperature/ 


500 K 


600 K 


700 K 


(Number of samples) 


(10 6 ) 


(10 6 ) 


(8.5xl0 5 ) 


Property 








£, 10~ 13 erg/atom .... 


-6.O04±.0Ol 


-5.602 ±.001 


-5.188 ±.001 


/M0 11 dyne/cm 2 


1.643 ±.001 


1.689 ±.001 


1.737 ±.001 


(3P/3T), 








10 7 dyne/cm 2 • K. . . 


5.0±.2 


4.8-1- .2 


4.7±.2 


C„, 10- 16 erg/K • atom. 


4.2 ±.2 


4.2 ± .2 


4.1 ±.2 


C n r , 10 n dyne/cm 2 


6.4 ±.3 


6.6 ±.3 


6.7±.3 


C 12 r , 10" dyne/cm 2 


0.97 ± .3 


1.1 ±.3 


1.1±.3 


C 44 r , 10" dyne/cm 2 


1.2±.2 


1.4 ±.2 


1.4±.2 


AC, 10 10 dyne/cm 2 


4.9 ±.5 


5.6 ±.5 


6.3±.5 
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FIGURE 2. Running estimates for the pressure (P), the derivative of the 
pressure with respect to temperature at constant volume ((3P/dT)y) and 
the isothermal elastic constant C, J for the T — 500 K state. 

Values are plotted for intervals of 5X10 4 trials. 

The estimates for the pressure and the energy are quite 
stable while those for the other quantities, all of which in- 
volve fluctuations, exhibit substantial variations with the 
number of Monte Carlo trials. The variation shown here is 
substantially larger than that found for the elastic constants 
of sodium by Cohen et al. [4]. 

4. Conclusions 

Our principal conclusion is that the temperature depend- 
ence of the elastic constants, for fixed volume, is a useful 
further test of pair potentials which satisfy low temperature 
tests. By examining the temperature dependence under 
fixed volume conditions, it is possible to sample the full 
range of the effective potential without introducing the 
complication of volume dependent changes in the potential 
itself. 

The calculated temperature dependence of the elastic 
constants is zero within the uncertainties quoted in table 2. 



This is certainly not the case for aluminum [17] where the 
elastic constants are decreasing functions of temperature. 
We see that the potential of figure 1, which works quite 
nicely at low temperatures is not at all a suitable potential 
function for aluminum. We also observe that a good deal of 
care is needed to insure that the Monte Carlo averages are 
stable. To achieve this, fairly long (in terms of the number 
of trials) simulations may be necessary. Finally, we conjec- 
ture that the local maximum shown in figure 1. between the 
first and second neighbor distances is responsible for the 
unphysical temperature dependence of the elastic con- 
stants. This conjecture is based on the idea that the max- 
imum inhibits anharmonic effects and thereby makes the 
lattice "stiffer" than it should be. 
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